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ABSTRACT 

The intrinsic distribution of spectral indices in GeV energies of gamma-ray-loud blazars is a critical 
input in determining the spectral shape of the unresolved blazar contribution to the diffuse extra- 
galactic gamma-ray background, as well as an important test of blazar emission theories. We present 
a maximum-likelihood method of determining the intrinsic spectral index distribution (ISID) of a 
population of 7-ray emitters which accounts for error in measurement of individual spectral indices, 
and we apply it to EGRET blazars. We find that the most likely Gaussian ISID for EGRET blazars 
has a mean of 2.27 and a standard deviation of 0.20. We additionally find some indication that FS- 
RQs and BL Lacs may have different ISIDs (with BL Lacs being harder). We also test for spectral 
index hardening associated with blazar variability for which we find no evidence. Finally, we produce 
simulated GLAST spectral index datasets and perform the same analyses. With improved statistics 
due to the much larger number of resolvable blazars, GLAST data will help us determine the ISIDs 
with much improved accuracy. Should any difference exist between the ISIDs of BL Lacs and FSRQs 
or between the ISIDs of blazars in the quiescent and flaring states, GLAST data will be adequate to 
separate these ISIDs at a significance better than 3<r. 

Subject headings: galaxies: active - gamma rays: observations - gamma rays: theory 



1. INTRODUCTION 

The Energetic Gamma-Ray Experiment Telescope 
(EGRET) aboard the Compton Gamma-Ray Observa- 
tory observed the gamma-ray sky in energies between 
30MeV and - 10 GeV between 1991 and 2000. The 
third (and last) catalog of EGRET point sources (Hart- 
mann et al. 1999) included 271 resolved objects of which 
93 were identified, either confidently or potentially, as 
blazars (gamma-ray loud active galactic nuclei). Thus, 
blazars constitute the class of 7-ray emitters with the 
largest number of identified members. 

The term "blazar" is used to refer collectively to BL 
Lac objects and 7— ray loud flat spectrum radio quasars 
(FSRQs). Blazars are believed to be active galactic nu- 
clei (AGNs) with the jet aligned with our line-of-sight 
(Blandford & Konigl 1979; see Urry & Padovani 1995 
for a review of general properties of blazars). The main 
contribution to the 7-ray emission of blazars is generally 
thought to result from cither inverse Compton scatter- 
ing by relativistic electrons in the jet of lower energy pho- 
tons, produced either by synchrotron emission within the 
jet (SSC), or by emission external to the jet, for example 
from the accretion disk around the central engine (EC); 
or from proton-induced cascades (see, e.g., von Montigny 
et al. 1995, Bottcher 2006 and references therein for a re- 
view of blazar 7— ray emission processes). 

Blazars fainter than the ones in the 3rd EGRET (3EG) 
catalog (and thus unresolved by EGRET) are expected 
to have a sizable contribution to the diffuse, isotropic 
gamma-ray background detected by EGRET (Sreekumar 
et al. 1998) which is presumably of extragalactic origin. 
The exact amount of this contribution remains unclear, 
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as different models for the blazar 7-ray luminosity func- 
tion result in very different predictions for the cumulative 
7-ray flux from all unresolved blazars, ranging from a few 
to 100% of the EGRB (Padovani et al. 1993; Stecker et 
al. 1993; Salamon & Stecker 1994; Chiang et al. 1995; 
Stecker & Salamon 1996a, hereafter SS96a; Kazanas & 
Perlman 1997; Chiang & Mukherjee 1998; Mukherjee & 
Chiang 1999; Miicke & Pohl 2000; Kneiske & Mannheim 
2005; Dermer 2007; Giommi et al. 2006; Narumoto & 
Totani 2006). Such uncertainties aside, blazars must be 
accounted for in any attempt to understand, model, or 
otherwise make use of the high-energy diffuse photon in- 
ventory (e.g. in constraining exotic high-energy physics). 
The lowest reasonable value we can expect for the cumu- 
lative diffuse intensity of unresolved faint blazars yields 
a strong lower limit for the level of the extragalactic 
gamma-ray background (EGRB). To constrain contribu- 
tions from any other plausible source, this lower limit 
must be subtracted from the observed background. 

A second critical property of blazars is their spectral 
index distribution (SID). The energy spectrum of the 7- 
ray emission of blazars in the EGRET energy range can 
be well approximated by a power law, F(E) oc E~ a with 
values of a, the spectral index, for individual blazars fit- 
ted to be in most cases between 2 and 3. Blazar spectra 
in the 7-ray regime encode important information con- 
cerning particle acceleration and emission processes in 
blazar jets (see von Montigny et al. 1995 and references 
therein). In addition, the distribution of blazar 7-ray 
spectral indices is a critical input in the estimation of 
the unresolved blazar contribution to the extragalactic 
diffuse background (SS96a; Pohl et al. 1997). Whether 
the collective emission of unresolved blazars may be the 
dominant component of the diffuse extragalactic gamma- 
ray background depends not only on intensity, but also 
on spectral shape (Stecker & Salamon 1996a, b; Pohl et 
al. 1997; Strong et al. 2004). Additionally, blazars may 
constitute the dominant component of the extragalactic 
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gamma-ray background only at certain energies (Pavli- 
dou & Fields 2002). In order to assess these issues, the 
spectral features of the unresolved blazar emission need 
to be studied, and hence the spectral index distribution 
of blazars needs to be understood. 

1.1. Past Work 

Obtaining the spectral index distribution of blazars 
presents three major difficulties: 

(1) Large measurement uncertainties in individual 
blazar spectral indices, due to low photon statistics. 
These errors contaminate the sampling of the underlying 
intrinsic spectral index distribution (ISID) 4 , by exagger- 
ating its spread, possibly to a large degree. 

(2) Possible systematic change of the spectral index 
with flaring [suggested, e.g., by von Montigny et al. 1995; 
Mukherjce ct al. 1996 (however, note that for the case 
of PKS 0528+134, the significance of the result was re- 
duced by subsequently obtained data, see Mukherjee et 
al. 1997); SS96a; Miicke et al. 1996; Pohl et al. 1997]. 

(3) Possible existence of two spectrally distinct popula- 
tions (BL Lacs and FSRQs) in the resolved blazar sample 
(e.g. Mukherjee et al. 1997; Pohl et al. 1997) 

Despite these difficulties, several authors have stud- 
ied different aspects of the statistical properties of GeV 
blazar spectra. SS96a calculated a spectral index dis- 
tribution for resolved EGRET blazars which they then 
used to derive the spectral shape of the collective emis- 
sion from unresolved blazars. They recognized that an 
appreciable spread in blazar spectral indices will lead to 
a pronounced concavity in the collective emission spec- 
trum (see Brecher & Burbidge 1972), a feature which is 
also at least tentatively present in determinations of the 
EGRB based on EGRET data (Sreekumar et al. 1998; 
Strong et al. 2004). Additionally, they insightfully em- 
phasized the potential importance of both variability as 
well as measurement errors in the determination of the 
blazar SID. SS96a remains to this day the only work as- 
sociating the distribution of blazar spectral indices as 
measured in individual objects with a model predicting 
the spectral shape of the unresolved blazar emission. 

However, their treatment suffers from three major un- 
resolved problems. First, their treatment of variability 
was based on very uncertain information from very few 
objects. Second, BL Lacs and FSRQs are treated as a 
single population, while the validity of this assumption 
was not evaluated. Finally, their treatment of measure- 
ment errors worsens, rather than alleviates, the overes- 
timation of the SID spread, and thus, overestimates the 
curvature of the unresolved blazar emission spectrum. 
This problem will not automatically disappear when the 
much larger sample of detected blazars and correspond- 
ing measured spectral indices of the upcoming GLAST 
mission becomes available. The bulk number of blazar 
detections will always be close to the instrument sensi- 
tivity limit, and will involve only a few tens of photons 
from each object, which means that the bulk number of 
measured spectral indices will always have substantial 
measurement errors associated with them. Therefore, 

4 In this paper, we will use the term "spectral index distribution" 
(SID) to refer to the distribution of measured spectral indices and 
the term "intrinsic spectral index distribution" (ISID) to refer to 
the true distribution of spectral indices of the blazar population 
(i.e. free of any contamination due to measurement errors). 



in order to be able to utilize all future measurements, 
including those for faint objects, in understanding the 
spectral properties of 7-ray-loud AGN in the GeV energy 
range, spectral index uncertainties have to be carefully 
dealt with. 

Pohl et al. (1997) used a different method to derive 
the spectral shape of the collective unresolved emission 
of blazars, which effectively circumvents the problem of 
large uncertainties in the measurement of spectral in- 
dices of individual objects. They co-added the spec- 
tra of resolved AGN and pointed out that if the spec- 
tral properties of unresolved blazars are similar to those 
of resolved blazars, then this co-added spectrum is the 
most appropriate quantity (in terms of spectral shape) 
for comparison with observations of the EGRB. They 
performed the analysis separately for BL Lacs and FS- 
RQs, and found that the BL Lac co-added spectrum is 
harder by 5a = 0.12 ± 0.08 than that of FSRQs, al- 
though the significance of the difference is low due to the 
low-number BL Lac statistics. Similarly, they found a 
difference between the FSRQ co-added flaring spectrum 
and the time-averaged co-added FSRQ spectrum, with 
the flaring spectrum being harder by 5a — 0.18 ± 0.05, 
a result with a statistical significance between 3 and 4c 
They verified that the co-added spectra are concave, al- 
though their method yields a much smaller spectrum cur- 
vature than the SS96a model. 

Mukherjee ct al. (1997), as part of a comprehensive 
analysis of the EGRET blazar observations available at 
that time, evaluated the average spectral properties of 
the blazar population and tested for indications of evo- 
lution of the spectral index with redshift and spectral dif- 
ferences between FSRQs and BL Lacs. Assuming equal 
errors in individual measurements of spectral indices, 
they found an average spectral index of 2. 15 ±0.04 for all 
blazars, and average spectral indices of 2.03 ± 0.09 and 
2.20±0.05 for the subsets of BL Lacs and FSRQs, respec- 
tively. They concluded that the statistical significance of 
the difference between the mean spectral indices of the 
two populations was at the level of 2. bo- and was not 
enough to justify a spectral separation of the two pop- 
ulations. Finally, they tested for redshift evolution of 
blazar spectra and found that for both BL Lacs and FS- 
RQs, the data were consistent with no evolution. Neither 
Pohl ct al. (1997) nor Muhkerjee et al. (1997) derived an 
SID for blazars. 

1.2. This Work 

The Gamma-ray Large Area Space Telescope 
(GLAST), which is scheduled for launch in 2007, 
will be able to address and ameliorate these difficulties. 
It will detect many more blazars than EGRET (between 
1,000 and 10,000 blazars; Stecker & Salamon 1999; 
Narumoto & Totani 2006; Dcrmcr 2007), enabling 
confirmation or rejection of any statistical trends seen 
in EGRET data. In order to facilitate the spectral 
studies of blazars with GLAST, it is important to use 
EGRET data to identify open questions and issues 
which can benefit from the dramatically improved 
GLAST statistics. Furthermore, it is advantageous 
to use EGRET data to identify and test appropriate 
statistical techniques which will allow us to deal with 
observational difficulties which are also expected to be 
present in GLAST data (such as large, varying errors of 
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measurement in individual spectral indices). 

In this context, with this paper we attempt to: (a) as- 
sess whether there is statistically significant evidence for 
evolutionary effects on the spectral index, manifesting as 
a correlation between the spectral index with cither red- 
shift or (isotropic) photon luminosity ($2]); (b) estimate 
the extent to which individual measurement errors affect 
the sampling of the SID of EGRET blazars ($3) and per- 
form a maximum likelihood analysis which accounts for 
these errors and determines the "most likely" parameters 
of the ISID (fj4j ; (c) assess whether there is statistically 
significant evidence for a difference between the ISIDs 
of BL Lacs and FSRQs ((JSJl and the ISIDs of flaring and 
quiescent blazars (©; (d) predict to what extent the up- 
coming GLAST observations will help us in resolving the 
issues mentioned above (jf7|). We summarize and discuss 
our conclusions in 

The dataset used for our analysis, except where explic- 
itly stated otherwise, is the set of 66 blazars characterized 
as "confident AGN identifications" in the 3EG catalog. 
We divide the population into 14 BL Lacs and 51 FSRQs 
(four sources have not been used due to lack of measured 
P1234 fluxes) as in Nolan et al. (2003). The spectral in- 
dices of individual objects are the average, P1234 indices 
quoted in 3EG. 

2. DOES THE BLAZAR SID DEPEND ON REDSHIFT OR 
LUMINOSITY? 

If blazars evolve spectrally with redshift, or if their 
spectral properties depend on luminosity, then the SID 
calculated from resolved, low-z and/or high-luminosity 
blazars is not representative of the spectral properties 
of unresolved, high-z and/or low-luminosity blazars, and 
should not be used to calculate the spectral shape of their 
collective emission. For this reason, we test for correla- 
tions between blazar spectral index and redshift or lumi- 
nosity. Figure [1] shows plots of the spectral index ver- 
sus photon luminosity (upper panel) and redshift (lower 
panel) for the 66 3EG confident blazars. The (isotropic) 
photon luminosity L p (shown in Fig. [1] in units of 10 50 
photons/sec) was calculated from the PI 234 photon flux 
Fpi234 using 



47rFp 1234 d|/(l + z). 



(1) 



Here dh is the luminosity distance, calculated from z 
assuming a concordance f2 m = 0.3, Q\ = 0.7, H n = 
70km s _1 Mpc _1 cosmology. 

Though there is no obvious trend in either panel of Fig. 
[T] we apply a formal non-parametric Spearman test for 
correlations (e.g. Wall & Jenkins 2003). Given a sample 
of N data pairs of variables, the two variables are ranked 
such that (Xi,Yi) are the ranks of the variables for the 
ith pair (1 < Xi < N and 1 < Yj, < N). Then, the 
Spearman rank correlation coefficient is computed: 



££i(*«-*) s 

N 3 - N 



= 1-6 



(2) 



with range < |r s | < 1; a high value indicates a 
significant correlation. The coefficient values we find 
for the spectral index/luminosity pairs are r s = 0.145 
(BL Lacs) and r s = —0.109 (FSRQs) and for the spec- 
tral index/redshift pairs are r s = 0.141 (BL Lacs) and 
r s = 0.238 (FSRQs). In all cases, such values of \r s \ 
or smaller occur by chance more than 10% of the time. 
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Fig. 1. — Spectral index versus luminosity (upper panel) and 
versus redshift (lower panel) for BL Lacs (filled circles/solid lines) 
and FSRQs (open circles/dashed lines). 



Therefore, we find no evidence for correlations between 
spectral index and either redshift or luminosity. Accord- 
ingly, we find no evidence for cosmological evolution in 
the particle acceleration and 7— ray emisison mechanisms 
operating in these AGNs, at least in the redshifts under 
consideration. 

3. THE EFFECT OF MEASUREMENT ERRORS 

We now turn our attention to the reconstruction of 
the ISID of blazars using measurements of the spectral 
indices for individual members of the blazar population. 
We assume that the ISID (probability that a randomly 
selected blazar has a true spectral index between a and 
a + da) can be adequately described by a Gaussian, 



Pintrinsic(a)rfa 



1 



V2vrcr 



exp 



(a - a ) 5 



2al 



da . (3) 



The spread of the ISID (which is equal to ctoj and should 
not be confused with the "error on the mean," the un- 
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certainty in our knowledge of ao) is not necessarily well- 
approximated by the spread of the distribution of mea- 
sured spectral indices. In sampling a distribution with 
significant errors in individual measurements, the spread 
of the resulting measured distribution may be dominated 
by the errors, especially if they are comparable with the 
width of the ISID. In this case, if the measured SID is 
used to calculate the spectral shape of the collective emis- 
sion of unresolved blazars, sources with spectral indices 
away from the mean would be overestimated and ulti- 
mately lead to an overestimate of the curvature of the 
collective emission spectrum. It is therefore necessary 
to determine the degree to which the measured SID is 
error-dominated before we can decide whether it is rep- 
resentative of the ISID. 

In the formalism implemented by SS96a, the measured 
SID is reconstructed by co-adding all of the probability 
density functions (assuming Gaussian errors) of individ- 
ual spectral index measurements of blazars in the second 
EGRET catalog: 



Pobs 



i(a)da 



1 

N 



N 

E 

»=i 



i 



2tt 



da, (4) 



where is the spectral index of blazar i, is the error in 
measurement of a,, and N is the total number of blazars. 5 
Equation ([4| is a good representation of the distribution 
of measured spectral indices (see Fig. [2] SS96a SIDs rep- 
resented by solid lines) and has the additional advantage 
that it does not make any a priori assumptions about the 
shape of the SID, but it could still be contaminated by 
measurement errors and as such, might not be a fair rep- 
resentation of the ISID. With this concern, we assess the 
significance of error in the SID, by performing a simple 
Monte Carlo test. 

We consider separately the samples of confident 
EGRET BL Lacs and FSRQs assuming that the ISID 
of each subset has the form of Eq. ([3]). We take the 
mean of each sample, ao, to be the "trimmed mean" 
(see e.g. Wilcox 1997) of the sample and the variance <7q 
to be variable. For later comparison, we also calculate, 
the "trimmed variance," of, of each population. We cal- 
culate ao and of by trimming the top and bottom five 
percent of the combined populations. In this way, we 
obtain for the BL Lacs: ao = 2.20 and at = 0.33; and 
for the FSRQs: a = 2.39 and a t = 0.22. 

We perform a number of "mock observations" from 
each ISID equal to the number of corresponding objects 
EGRET detected, with a randomly chosen 3EG uncer- 
tainty for each object, and calculate the trimmed vari- 
ances of each set. If some assumed ISID spread, sampled 
with EGRET uncertainties, more frequently results in 
simulated trimmed variances smaller (larger) than that 
of the corresponding EGRET dataset, then it is most 
likely too small (large) compared with the ISID spread 
occurring in nature. The spread of the ISIDs of BL Lacs 
and FSRQs occurring in nature should be comparable 
to the spread for which the simulated datasets most fre- 
quently have trimmed variances comparable to those of 

5 SS96a also accounted for the different blazar states (flar- 
ing/quiescent) by shifting this distribution (toward harder/softer 
indices), but they considered BL Lacs and FSRQs as a single pop- 
ulation. 
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Fig. 2. — Top: Histogram of spectral indices of 3EG BL Lacs 
with SS SID (solid) overlayed. Bottom: Same for the 3EG FSRQs. 
The dashed lines represent the maximum-likelihood Gaussians for 
the same sets (see ifll . 



the EGRET datasets. Figure [3] demonstrates this com- 
parison between the trimmed variances of the EGRET 
set and the simulated sets. The fraction of simulated 
sets with trimmed variance greater than that of the cor- 
responding EGRET set is plotted against the spread of 
the parent ISID. The median ISID spread is equal to 
0.27 for the BL Lacs and 0.20 for FSRQs. If we were to 
fit a Gaussian to the SS96a prescription for the SID, we 
would obtain a standard deviation of 0.46 in the case of 
BL Lacs and 0.36 in the case of FSRQs - almost twice the 
preferred values of the Monte-Carlo analysis. Hence, we 
conclude that the SS96a prescription does overestimate 
significantly the spread of the ISID. 



THE ISID OF EGRET BLAZARS 
APPROACH 



A LIKELIHOOD 



Having shown that measurement errors can signifi- 
cantly contaminate the determination of the ISID, we 
now employ a likelihood analysis which allows us to ex- 
plicitly account for these errors and constrain the pa- 
rameters of the ISID. Explicitly accounting for measure- 
ment errors is necessary as ignoring them may not sim- 
ply increase the uncertainty in our estimated parameters 
but lead to incorrect parameter inferences (e.g. Loredo 
2004). 

Given a set of parameters Xi , which define a statistical 
distribution, and a dataset yj , the probability of Xi hav- 
ing certain values given the data is proportional to the 
probability of measuring yj given those values of Xi (the 
likelihood) times the probability that Xi has those partic- 
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Fig. 3. — Fraction of sets for which the trimmed variance of 
the set with errors is greater than that of the EGRET set without 
errors for sets with the same number of objects as confident BL 
Lacs (solid line) or FSRQs (dotted line). The dot-dashed lines 
indicate the median "true spread" of the simulated samples, which 
is equal to 0.27 for BL Lacs and 0.2 for FSRQs. 



ular values before measurements are taken (the prior): 



P(xi\yj) oc P(xi) x C(yj\xi). 



(5) 



If the prior is flat (constant for all values of Xi), maximiz- 
ing the likelihood allows us to determine the most prob- 
able values for the parameters Xi (see, e.g., Lee 1989). 

In our analysis, we assume that both the ISID and the 
errors are Gaussian. The parameters we wish to calculate 
are the mean spectral index, ao, and spread, ao- Without 
error, the likelihood of measuring a spectral index a given 
a Gaussian ISID with a mean of ao and a spread of ao 
would be I — exp [—(a — ao) 2 /2a 2 ] /V2~jrao- To include 
measurement error, we need to distinguish between the 
true spectral index of an object, a (which is unknown and 
is therefore a nuisance parameter over which we need to 
marginalize), and its measured spectral index, ctj. Then, 
the likelihood for a single spectral index measurement 
becomes 



da 



exp [-(a - a,) 2 /(2a])] exp [-(a - a ) 2 /(2a 2 )] 



(6) 



2nai 



2nao 



where the subscript denotes blazar j. If we have N such 
independent spectral index measurements, then the over- 



all likelihood becomes L = J| 
normalization factors, 



JY 



= 1 3 



lj or, omitting constant 



C 



N 

n 



1 



a 2 + a 2 



exp 



1 N 



(a 3 -a ) 2 



(7) 



The maximum-likelihood mean and spread of the ISID 
can be found by maximizing the likelihood; hence, by 
simultaneously solving the equations (for derivation, see 
Appendix) 



a 



E^S E^- 

^ a 2 + a 2 ^ a 2 + a 



(8) 
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Fig. 4. — Likelihood function contours for the Mattox 2001 (solid 
lines) and the 3EG (dashed lines) confident blazar samples. The 
l<x, 2tx, and 3(T contours we plot throughout this paper represent 
equal-likelihood contours which include 68%, 95.5% and 99.7% of 
the total volume under the likelihood surface. The x and y axes 
represent the mean (»o) an d spread (<tq) of the ISID respectively. 
The maximum-likelihood parameters of each ISID are denoted by 
x. 



and 



N 



N 

E 



(qj - ap) 2 

.2X2 



(9) 



In Fig. 01 we plot likelihood contours for the confident 
blazar set of Mattox 2001 (solid lines) and of the 3EG 
(dashed lines). The maximum-likelihood parameters for 
the ISID in the case of the Mattox set are ao — 2.27 and 
ao = 0.20, and those for the 3EG set are ao = 2.29 and 
ao = 0.22. As we can see in Fig. 21 the likelihood func- 
tions of the two sets are consistent with each other, which 
shows that our analysis is not sensitive to the inclusion 
or exclusion of a few members. 

5. BL LACS VS FSRQS: DO THEIR ISIDS DIFFER? 

Gamma-ray-loud BL Lacs and FSRQs have different 
properties in the GeV energy range (e.g. different vari- 
ability properties, Vercellone et al 2004; different mean 
spectral properties, Mukherjee et al 1997, Pohl et al 1997; 
different redshift distributions and possibly different lu- 
minosity functions, Miicke & Pohl 2000, Dermer 2007). 
It is therefore reasonable to test whether the maximum- 
likelihood ISID for the BL Lacs and FSRQs differ at a 
statistically significant level. 

We apply the analysis of Sj4]to the sets of BL Lacs and 
FSRQs of the 3EG catalog. Fig. [5] shows the likelihood 
function contours for these populations. The ISID pa- 
rameters are found to be ao = 2.15, ao = 0.28 for BL 
Lacs and a = 2.3, a a = 0.19 for FSRQs. The likelihood 
functions of the two populations do indicate a marginal, 
la separation. This possible spectral differentiation be- 
tween the two populations, if true, will be confidently 
confirmed by GLAST data (see SJTJ) . The maximum- 
likelihood ao's are close to the trimmed means of each 
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Fig. 5. — Likelihood function contours (as in Fig. [4j for 3EG BL 
Lacs (solid lines) and FSRQs (dotted lines). 

dataset, while the maximum-likelihood oo's are almost 
identical to the median "true spread" for each set, as 
determined by the Monte Carlo test outlined in <J3] indi- 
cating the consistency in the results of the two analyses. 
Our maximum-likelihood Gaussians for BL Lacs and FS- 
RQs are overplotted on Fig. [2] with the dashed lines. 

6. BLAZAR VARIABILITY AND THE MEASURED 
SPECTRAL INDEX 

Though a blazar may spend most of its time in a quies- 
cent state, it can undergo periods of flaring during which 
its flux is significantly enhanced, up to an order of magni- 
tude (e.g. McLaughlin et al. 1996; Mukherjee et al. 1997; 
Nolan et al. 2003; Vercellone et al. 2004). This increase in 
flux may introduce a detection selection effect for fainter 
blazars favoring the flaring state. On the other hand, 
blazars are more likely to be pointed at during quies- 
cence since they spend more time in the quiescent state. 
The spectral index of a blazar is determined using the in- 
tegrated EGRET maps which can include photons from 
both states, thus complicating the determination of the 
ISID: it is unclear whether the spectral index distribution 
of blazars determined from EGRET data is representa- 
tive of blazars in their flaring or quiescent state. This 
distinction is important because it has been suggested 
that the blazar spectrum might change during the flaring 
state. If such a trend is indeed present and most photons 
involved in EGRET blazar detections come from blazars 
in the flaring state, then the derivation of the SID will be 
biased towards spectral indices representative of flaring 
blazars (as was assumed in SS96a). 

Determining whether a blazar is in a flaring or quies- 
cent state during a certain viewing period is made com- 
plicated by low photon statistics. During the long ex- 
posure necessary to detect a blazar, the source could be 
in both states, and the duty cycle is largely unknown 
(Vercellone et al. 2004). Furthermore, the source identi- 
fication and flux estimation for different viewing periods 
and for the cumulative map (from which the spectral in- 
dex is derived) is done by separate statistical processing 



of each map, and as a result, the photon counts from the 
individual viewing periods don't add up to those of the 
cumulative map. With these uncertainties in mind, we 
approach the problem using a simple recipe. 

In a manner similar to that of Vercellone et al. (2004), 
we classify each individual viewing period of a confi- 
dent AGN (excluding periods for which only flux upper 
limits were quoted) as a "mostly" flaring period if the 
flux of the blazar in the particular viewing period, Fyp, 
is > T x FP1234 (where T > 1 is some constant en- 
hancement factor) and if the error on the measurement 
of Fyp is less than the difference between the Fyp and 
-Ppi234; otherwise, we classify it as a "mostly" quiescent 
period. As in Vercellone et al. (2004), and in Jorstad et 
al. (2001), we consider the conventional case 6 of T = 1.5. 
We also consider the case of T = 3.75, chosen as the en- 
hancement factor for which ~ 75% of the objects have a 
flaring fraction below 0.2. 

Once we have classified all of the viewing periods dur- 
ing which a blazar was detected, we add all of the photon 
counts from the flaring periods to get the number of flar- 
ing photons, Nf, and likewise for the quiescent periods to 
get the number of quiescent photons, N q . We then cal- 
culate the "flaring photon fraction," / = Nf/(Nf + N q ), 
for every confident blazar. 

Figure [6] shows the fraction of blazars with a flaring 
fraction greater than / as a function of /, for both en- 
hancement factors. As expected, the flaring photon frac- 
tion depends on our definition of a "flaring state:" the 
fraction of objects with flaring photon fraction > /, is 
systematically lower for all values of / for the higher 
value of J- . However, even for T — 1.5, the median flar- 
ing photon fraction is not very high (equal to 0.5 for BL 
Lacs, and 0.6 for FSRQs, while only 20% of BL Lacs 
and 40% of FSRQs have a flaring photon fraction higher 
than 0.7, at which point photons predominantly come 
from a flaring state). Therefore, the photons that are 
used to derive individual blazar spectra represent a bal- 
anced mix of photons originating in flaring and quiescent 
states, if not photons that come primarily from the quies- 
cent state (with the quantitative details of this statement 
depending on the value of the enhancement factor). This 
analysis suggests that even if spectral indices of blazars 
do harden considerably during the flaring states for most 
objects, there is no indication that the photon budget is 
overly biased towards one variability state. Hence, if a 
SID is determined from EGRET data, there is no evi- 
dence that it would be predominantly representative of 
the flaring-state or quiescent-state SID. 

Qualitatively, this result also has implications for the 
blazar duty cycle. If the time spent flaring is compa- 
rable to the time spent in quiescence, we would expect 
the number of flaring photons to be significantly greater 
than the number of quiescent photons, yielding a flaring 
fraction close to one. This however is not the case even 

6 Note that Vercellone et al. (2004) and Jorstad et al. (2001) com- 
pared the individual viewing periods with the inverse-uncertainty- 
weighted mean of fluxes and upper limits of the individual view- 
ing periods. However, this weighted flux is more biased towards 
higher fluxes than the true long-term average flux of the blazar, 
as a blazar is preferentially detected in a viewing period when it 
is flaring (Vercellone et al. 2004). Thus, we systematically identify 
more "flaring" periods than Vercellone et al. (2004) and Jorstad et 
al. (2001). 
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Fig. 6. — Fraction of objects with fraction of flaring photons 
greater than / as a function of / for T = 1.5 (upper panel) and 
T = 3.75 (lower panel) for BL Lacs (solid lines) and FSRQs (dotted 
lines) . 



for a lower value of T . Therefore, these blazars must 
have been viewed more often in their quiescent state. 

6.1. Is there statistical evidence for spectral index 
hardening during flaring? 

As a first indication for spectral index shift with flaring, 
we can test for a spectral index/flaring photon fraction 
correlation. No obvious trend exists but, as in fj2) we 
also perform a non-parametric Spearman rank coefficient 
test, treating the BL Lacs and FSRQs separately, and for 
T = 1.5 and T — 3.75. 

With the exception of BL Lacs in the T = 3.75 case, 
the data are consistent with uncorrelated variables at the 
20% level. In the case of T = 3.75, all BL Lacs but one 
have / = 0, which gives a false positive result in the cor- 
relation test. However, when we treat the BL Lacs and 
the FSRQs as a single population, the results are again 
consistent with uncorrelated variables. A recent sys- 
tematic reanalysis of EGRET data for all blazars bright 
enough for time-resolved spectroscopy by Nandikotkur et 
al. (2007) has also yielded similar results despite using 
a radically different approach, as no systematic trend of 
the spectral index with changing flux was found. 

A second method to look for a systematic spectral in- 
dex shift with changing flux is to perform the maximum- 
likelihood analysis of <g] separately for the populations 
of "mostly flaring" and "mostly quiescent" blazars. The 
set of "mostly flaring" blazars are the 22 EGRET blazars 
that have a flaring photon fraction / > 0.7, for T = 1.5. 
The set of "mostly quiescent" blazars are the 10 EGRET 
blazars that have flaring photon fractions / < 0.3. Likeli- 
hood function contours for the ISID parameters of these 
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Fig. 7. — Likelihood function contours (as in Fig. [4]l for 
mostly flaring (solid lines) and mostly quiescent (dashed lines) 
blazars. Upper panel: 3EG blazars; lower panel: simulated 
EGRET dataset. 



two sets are shown in the upper panel of Fig. [71 The 
ISIDs of the two sets are fully consistent with each other, 
with maximum-likelihood parameters for the "mostly 
flaring" blazars a = 2.25 and cr = 0.18, and for the 
"mostly quiescent" blazars «o = 2.13 and Co = 0.11. If 
anything, the maximum-likelihood ISID of mostly flar- 
ing objects peaks at slightly softer indices than that of 
the mostly quiescent sources (a trend that could be theo- 
retically accomodated if, for example, softer EC emission 
overtakes an otherwise harder, SSC dominated, spectrum 
during flaring; see e.g. Bottcher 1999). 

In order to determine whether the EGRET data would 
be sufficient to reveal spectral index hardening with flar- 
ing if such a trend does exist, we perform the follow- 
ing test. We assume that blazar spectral indices do in- 
deed harden with flaring, and that their ISIDs differ as 
suggested by SS96a: they are identical Gaussians with 



Co = 0.20 and means that are determined by shifting the 
mean of the combined population, aa = 2.27 (mean and 
spread obtained from the maximum-likelihood analysis 
of the Mattox et al. 2001 set), by -0.05 for the flaring 
blazars and by +0.2 for quiescent blazars. We then sam- 
ple this dataset with EGRET uncertainties (as in S}3|) and 
create a mock EGRET dataset. We apply the same like- 
lihood analysis to these simulated datasets. Likelihood 
function contours are plotted in the lower panel of Fig. 
[7J Although the maximum-likelihood analysis correctly 
identifies the flaring population ISID as being harder 
than the quiescent population ISID, the separation of 
the maxima of the likelihood functions has a marginal 
la significance. Additionally, in the EGRET dataset the 
measurements of the spectral index do not come from 
pure "flaring" or "quiescent" states, and even the defini- 
tion of a "flaring" or "quiescent" state is dependent on 
the selected value of the enhancement factor T . Thus, 
we conclude that although EGRET data do not seem to 
indicate a hardening of the spectral index with flaring, 
they are not sufficient to exclude such a possibility either. 
However, as we will see in the next section, GLAST data 
will settle the issue. 

7. PROSPECTS FOR GLAST 

With the launch of GLAST later this year, the num- 
ber statistics of resolved blazars are expected to im- 
prove dramatically. Depending on models used to de- 
termine blazar numbers at lower fluxes, GLAST will de- 
tect between 1,000 - 10,000 blazars with flux > 2 x 
10~ 9 photons s _1 . However, most of these blazars will 
have fluxes close to the lower end of GLAST sensitivity 
with individual photon statistics and spectral index un- 
certainties comparable to those of the EGRET dataset. 
Thus, the extent to which the GLAST dataset will im- 
prove our understanding of the issues discussed in this 
paper is not immediately obvious. 

To assess this question, we construct simulated 
datasets of spectral indices and spectral index uncer- 
tainties for GLAST-detectable blazars. Based on these 
datasets, we predict by how much GLAST observations 
will improve the determination of statistical properties 
of the spectral indices of gamma-ray-loud blazars. For 
our study, we use the Dermer (2007) gamma-ray lumi- 
nosity functions for BL Lacs and FSRQs. These lumi- 
nosity functions are more conservative than those of, 
e.g., Narumoto & Totani (2006), Chiang & Mukherjee 
(2001) and SS96a as they predict a smaller number of 
detectable blazars; additionally, they treat BL Lacs and 
FSRQs separately, enabling us to distinguish between the 
two populations in our analysis. If GLAST detects more 
blazars than predicted in Dermer (2007) and assumed 
here, the statistics will only improve, strengthening the 
significance of our results. 

In assigning measurement uncertainties to the spectral 
index of each blazar in our simulated datasets, we make 
use of the strong anticorrelation of spectral index uncer- 
tainties in the EGRET dataset with the total number 
of detected P1234 photons. A single power-law can be 
fitted to the 3EG data, 

We expect that, despite any differences between the 
Large Area Telescope (LAT) and EGRET instruments, 



0.4 



0.3 : 




0.1 : : 

0.0 [ i i i : 

2.0 2.1 2.2 2.3 2.4 
«o 

Fig. 8. — Likelihood function contours (as in Fig. [4} for simulated 
GLAST datasets of BL Lacs (solid lines) and FSRQs (dashed lines) 

this empirical relation can also give a rough approxima- 
tion of the spectral index uncertainty for blazars detected 
by GLAST. 

7.1. Spectral separation between BL Lacs and FSRQs 

Using simulated GLAST datasets, we examine the pos- 
sibility that, if BL Lacs and FSRQs are indeed spectrally 
distinct populations, GLAST observations will be suffi- 
cient to separate them with a significance greater than 
3a. We assume that the ISIDs for BL Lacs and FSRQs 
are the maximum-likelihood Gaussians determined from 
EGRET data (a = 2.15 and a = 0.28 for BL Lacs and 
a = 2.3 and cr = 0.19 for FSRQs). We generate the 
simulated BL Lac and FSRQ datasets by appropriately 
sampling these ISIDs with GLAST uncertainties. The 
number of objects in each set and their fluxes are de- 
rived from Fig. 2 of Dermer (2007). From the flux of each 
object, we derive an approximate number of detectable 
photons, Aphoton,i = FiAtA, where Fi is the flux of the 
particular blazar, At is the estimated exposure time (as- 
suming that any given objects will be in the GLAST 
field of view for ~ 20% of a 2 year campaign), and A is 
the effective collecting area of the detector (8, 000cm 2 for 
GLAST). The spectral index uncertainty for each object 
is then derived from Eq. 1101 

Performing the analysis outlined in $4] on these 
datasets, we get the likelihood function contours shown 
in Fig. [5J The maximum-likelihood parameters of the 
ISIDs are «o = 2.17 and ao — 0.27 for the BL Lacs, 
and ao = 2.31 and ao = 0.2 for the FSRQs, excellent 
approximations of the assumed underlying ISID param- 
eters. The separation between the two populations is 
clearly established with significance greater than 3a. 

7.2. Spectral separation between flaring and quiescent 

blazars 

Finally, we investigate the possibility that GLAST ob- 
servations will be sufficient to determine whether flaring 
and quiescent blazars are indeed spectrally distinct. We 
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construct simulated GLAST datasets to test two scenar- 
ios: that the ISIDs of flaring and quiescent blazars have 
the same shape, but flaring blazars are shifted towards 
harder indices by 0.25 (as in $6.1$ ; and that flaring and 
quiescent blazars have ISIDs with parameters equal to 
the maximum-likelihood parameters yielded by the anal- 
ysis of the true EGRET dataset (note that in this case 
the flaring blazars are somewhat softer than quiescent 
blazars). 

As there are many more blazars in this analysis, we 
can account for the disparities between the flaring frac- 
tions for bright and faint blazars by calculating (from 
EGRET data) the fractions of bright and faint objects 
that are mostly flaring (/ > 0.7) or mostly quiescent 
(/ < 0.3); thereby, more carefully determining the num- 
bers of mostly flaring and mostly quiescent blazars. We 
assign true spectral indices as in jj6.II and spectral index 
uncertainties based on the flux as in £17.11 treating BL 
Lacs and FSRQs as a single population. 

We then apply the same likelihood analysis to these 
simulated datasets. Their likelihood function contours 
are plotted in Fig. [9] Although the expected partial con- 
tamination of any single state with both flaring and qui- 
escent photons may weaken this result, likelihood func- 
tions for the flaring and the quiescent blazars are well 
separated with a significance greater than 3a in both 
scenarios, strongly suggesting that the GLAST data will 
reveal any systematic change of spectral index with flar- 
ing, if present. 

8. DISCUSSION AND SUMMARY 

In this work, we have extensively studied possible 
trends in the spectral behavior of EGRET blazars, and 
have attempted to reconstruct the blazar ISID from 
EGRET observations, explicitly accounting for the large 
measurement errors arising from low-photon statistics. 
Additionally, we have identified still open questions on 
which progress is expected to be made once GLAST ob- 
servations become available, and we have estimated the 
extent of the expected improvement. Our conclusions 
can be summarized in the following statements: 

(a) We have investigated the possibility that source 
evolution could interfere with the application of our 
ISIDs to unresolved blazars. We have found no evidence 
for correlations between spectral index and redshift or 
between spectral index and luminosity in EGRET data. 
Therefore, we have concluded that there are no indica- 
tions for spectral index evolution in EGRET data. 

(b) We have used a Monte Carlo test to study 
whether measurement errors in individual spectral in- 
dices severely contaminate the observed SID. We found 
that measurement errors can indeed have a profound ef- 
fect and need to be accounted for properly. This will also 
be the case for GLAST data, as most of the newly de- 
tected blazars will suffer from similar low-photon statis- 
tics as with the EGRET blazars. We have performed a 
likelihood analysis which explicitly accounts for signifi- 
cant and unequal errors of measurement, and we found a 
most likely Gaussian ISID for the population of EGRET 
blazars with a mean of 2.27 and a spread of 0.20, signif- 
icantly narrower than the observed SID. 

(c) We have derived separate ISIDs for BL Lacs and 
FSRQs finding that FSRQs are well represented by a 
Gaussian with mean of 2.3 and spread of 0.19, while the 
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Fig. 9. — Likelihood function contours as in Fig. [7]for simulated 
GLAST data. Upper panel: likelihood contours assuming flaring 
blazars are harder; lower panel: likelihood contours assuming flar- 
ing and quiescent blazars ISID with parameters determined in the 
analysis of the EGRET simulated variability dataset. 



BL Lacs are better represented by a Gaussian with mean 
of 2.15 and spread of 0.28. Based on EGRET data, the 
significance of the spectral separation of the two popula- 
tions is marginal. We determined whether each spectral 
index of the EGRET confident blazars was more rep- 
resentative of the flaring or quiescent state by calculat- 
ing the flaring photon fraction for each object. We per- 
formed statistical tests for spectral index hardening and 
found no evidence for this effect. We calculated separate 
maximum-likelihood ISIDs for mostly flaring and mostly 
quiescent objects finding no evidence for a systematic 
separation of the two populations. However, we also 
showed that such a separation can neither be excluded 
based on EGRET data. This result is consistent with the 
recent findings of Nandikotkur et al. (2007) who stud- 
ied spectral index shifting during flaring on an object- 
to-object basis and also found no systematic trend. Ul- 
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timately, the question of hardening will be addressed by 
GLAST through more confident measurements of spec- 
tral indices and an increased ability for time-resolved 
spectroscopy (Dermer & Dingus 2004) allowing confident 
measurements of spectral indices for separate variability 
states of blazars. 

(d) We have made explicit estimates about the im- 
provement of our statistical analysis once GLAST data 
become available finding that, with the detection of many 
more blazars, GLAST will be able to separate with sig- 
nificance > 3<7 the ISIDs of BL Lacs and FSRQs, if these 
populations are indeed spectrally distinct. Similarly, we 
have examined two different scenarios for a possible sys- 
tematic spectral shift during flaring finding that in both 
cases, GLAST would be able to distinguish between the 
two ISIDs at > 3<T. 

The maximum-likelihood analysis presented in this 
work is therefore a useful tool in maximizing the infor- 
mation that can be obtained from the high-uncertainty 
measurements of spectral indices in the GeV band, where 
photon statistics for most objects are very low. Addition- 



ally, it can provide estimates for both the mean and the 
spread of the underlying, intrinsic distribution(s) of spec- 
tral indices of 7-ray emitters. Both of these aspects of 
the ISID(s) are not only important tests for 7-ray emis- 
sion models for blazars, but also essential inputs in de- 
termining the spectral shape of the diffuse emission from 
unresolved blazars (Pavlidou & Venters 2007). 
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Universite de Paris 7 for their hospitality and cama- 
raderie during the final stages of this work. This work 
was supported by the Kavli Institute for Cosmological 
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APPENDIX 



DERIVATION AND MAXIMIZATION OF LIKELIHOOD FUNCTION FOR THE ISID PARAMETERS 

As discussed in SJ31 the likelihood of observing N spectral indices aj (j = 1,...,N) with individual measurement 
uncertainties aj (j = 1, .., N) (where the aj are assumed Gaussian) if the ISID is Gaussian with mean ckq and variance 

(Tq is C = n^=i hi w ith 
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we can rewrite the exponentials in Eq. (IA1|) by completing the square as 
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Substituting this in Eq. (IA1|) and performing the integration we obtain 
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The likelihood function then becomes 
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Local maxima of this function will satisfy 
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Equation (|A5[) implies 
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therefore the first of Eqs. (|A6j) becomes 
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where we have used the identity 
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which is straight-forward to prove by induction. Therefore, the second of Eqs. (|A6[) becomes 
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